Global minimum potential energy conformations of small molecules

1. I n t r o d u c t i o n The s t u d y of molecu la r confo rmat ions is as fasc ina t ing a sub jec t as i t is u t t e r ly complex . A l t h o u g h the bas ic bu i ld ing blocks of molecules , the a toms , r e m a i n virt ua l l y unchanged in different c o m p o u n d s the versa t i l i ty of the ways t ha t they can be combined and reconfigure the o b t a i n e d a t o m chains resul ts in m a n y different spa t i a l conf igura t ions for a given molecule . One of these conf igurat ions , the mos t s tab le one, is of p a r t i c u l a r i m p o r t a n c e because i t d ic ta tes mos t of the p rope r t i e s of the molecule . Th i s p rovides enough m o t i v a t i o n to pred ic t the mos t s tab le confo rma t ion of a molecule based solely on the energet ics of the in te rac t ions be tween the a t o m s compos ing the molecule . Molecular mechanics is a widely used m e t h o d designed to provide a pr ior i accura te r epresen ta t ions of s t ruc tures and energies for molecules . I t o r ig ina t ed in 1946 when Hill [13] p roposed t ha t van der W a a l s in te rac t ions a long wi th s t r e tch ing and bend ing de fo rma t ions can be used for express ing the po t en t i a l energy of a molecule . Dos t rovsky , Hughes and Ingold [9, 7] u t i l ized this s ame bas ic pr inc ip le in an effort to u n d e r s t a n d the ra tes a t which var ious ha l ides underwent the S N 2 reac t ion . However, i t was Wes t ehe imer and Mayer [34, 24] who first used mo lecu l a r mechanics ideas successfully in exp la in ing the r a t e of r a c e miz a t i on of *AUTHOR TO WHOM ALL CORRESPONDENCE SHOULD BE ADDRESSED. 136 M A R . A N A S A N D F L O U D A S some optically active compounds. With the advent of computers during the 1950% molecular mechanics became one of the standard methods of structural chemistry. The most accurate representation of the potential energy of a molecule is the ab initio quantum mechanical approach. Quantum mechanical calculations, however, for all but the simplest molecules such as H2 or HCl are not feasible because of the large associated computational effort. As a result, tractable potential energy expressions had to be derived which adequately captured the energy contributions resulting from various types of atom interactions. In doing so, trade-offs between maximizing the accuracy of the potential energy expression and at the same time minimizing the computational effort to evaluate these functions had to be addressed. In quantum mechanics the Born-Oppenheimer approximation is routinely used. It states that the SchrSdinger equation for the molecule can be separated into a part describing the motions of the electrons and a part that describes the motion of the nuclei; and thus these two sets of motions can be studied independently. This means that the potential energy of a molecule in the ground electronic state is a function of only the nuclear positions. This defines a potential energy multi-dimensionM hypersurface in the coordinate set of the nuclear positions. In molecular mechanics this hypersurface is simply called potential energy surface and encompasses the effect on the potential energy of all possible conformations that a molecule can resume. Any point of this surface corresponds to a different conformation of the molecule and local minimum points on this surface are referred to as conformers. The systematic identification of the global minimum point on this surface is the focus of this work. Experimental evidence [28] shows that in the great majori ty of cases this total potential energy global minimum point corresponds to the most stable conformation of the molecule. However, there exists examples where this is not true due to the interplay of rotational-vibrational motions [25]. Molecular mechanics calculations employ an empirically derived set of potential energy contributions for approximating the Born-Oppenheimer surface [3]. This set of potential energy contributions, called the force field, contains adjustable parameters that are selected in a such a way as to provide the best possible agreement with experimental data. The main assumption introduced in molecular mechanics is that every parameter is associated with a specific interaction rather than a specific molecule. These parameters can be bond lengths; covalent bond angles; bond stretching, bending, or rotating constants; non-bonded atom interaction constants, etc. Thus, whenever a specific interaction is present, the same value for the parameter can be used even if this interaction occurs in different molecules [14]. Although it is not possible to prove the validity of this assumption, experimental results provide sufficient evidence that it is a reasonable assumption in most cases. Molecular mechanics force fields include two-body, three-body or even n -body interaction terms. v : E + E + E + . . . i>j i>j>k i>j>k>l It can be seen that the total potential energy V is not simply the sum of all pair interactions; it may involve energy terms which depend on the position and propGLOBAL MINIMUM ENERGY MOLECULAR CONFORMATIONS 137 erties of three or even more atoms due to the fact that the charges on the particles are in general polarized. However, for nonpolar atoms these terms can be neglected with little effect on the total potential energy [14]. This means that for nonpolar molecules the total potential energy expression is given approximately by, V = ~ V ~ j ( B ) + ~--~V~j(NB) i > j i > j where B indicates bonded atom interactions and NB nonbonded atom interactions. In most conformational calculations the bonded atom interactions are assumed to be constant and independent of the actual conformation and thus they can be omitted since they do not contribute to the configuration of the molecule. Also, in general only pairwise additive terms are significant in conformational calculations although higher order terms can become important under certain conditions. The general behavior of two-body potentials is very well studied. The force field is repulsive due to coulombic nuclear-nuclear and electron-electron interactions and attractive due to electron-nuclear interactions. In practice the pairwise potential interactions are usually represented by Lennard-Jones, Buckingham's, or Kihara's classical empirical potential functions: 1. L e n n a r d J o n e s 6 1 2 P o t e n t i a l F u n c t i o n Bij Aij ~)ij ~ r l 2 r6 The attractive term A i j / r 6 is rigorously derived for a pair of identical spherically symmetrical and chemically saturated molecules, and Bij /r 12 is an approximation for the repulsive component for which no rigorously derived expression exists. 2. B u c k i n g h a m P o t e n t i a l F u n c t i o n c~j vij = Aij exp (Bijrij) + r~j The Buckingham potential function differs from the Lennard-Jones function only in the form of the repulsive term. A two-parameter exponential function is used which in principle should be more specific in describing the repulsive interaction because of the additional parameter. K i h a r a P o t e n t i a l F u n c t i o n V i j ---4Co r i J / ~ 7 r i j / ~ With the help of the parameters 7, a Kihara's potential introduces an effective core and a shape dependence. 3. 138 M A R A N A S A N D F L O U D A S Before solving the resulting energy minimization problem it is necessary to decide whether or not the energy minimization will be performed on an independent set of internal coordinates (all bond lengths, covalent bond angles, and torsion angles) or on the Cartesian coordinates. Because it is easier to calculate the internal coordinates which are needed for evaluating the potential function from Cartesian coordinates rather than from an independent set of internal coordinates, most minimization methods use Cartesian coordinates. A number of methods have been proposed for finding the most stable conformation of a molecule through the identification of the global minimum point of the potential energy surface. All methods a t tempt to locate this point by tracing paths on the potential energy surface conjecturing that some of them will converge to the global minimum point. There are two main groups of methods; simulation type methods and gradient methods. Simulation type methods, including Monte-Carlo minimization and simulated annealing share the problem of selecting good strategies for "temperature" reduction dependence, optimal step size selection, and efficient random generator algorithms for generating the random walk. Gradient type methods are the most widespread methods for potential energy minimization of molecules. They can be divided into first derivative and second derivative minimization techniques. First derivative techniques (steepest descent) [35] follow a path defined by the steepest descent direction at every point. They perform satisfactorily only if all first derivatives are of the same order of magnitude, otherwise scaling problems cause oscillations. It has been also reported [2] that torsion angles are often not well minimized with these methods. Second derivative techniques, although more complicated, are far superior to first order methods. By utilizing second order derivative information calculated either analytically [1, 33, 29] or numerically [16], improved convergence rates are achieved. It is important to note that most energy minimization procedures do not locate energy minima, but rather stationary points which occur when all first order derivatives are equal to zero. Therefore such procedures may converge to a saddlepoint or even a energy maximum on the potential energy surface rather than to an energy minimum. For example, in the case of cyclohexane the boat conformation (C2,) or the conformation that has all six carbon atoms on a plane (D6h) may erroneously appear as energy minima. However, the main limitation shared by all the aforementioned methods is that unless there is a single potential well, the obtained minimum energy conformation depends heavily on the supplied initial conformation. This is not surprising because all currently available methods are local optimization methods guaranteed to find a local minimum at best. This is why in practice many trial geometries serve as initial points for the employed optimization method. These geometries are usually chosen from Dreiding models, or other similar considerations and thus there is no guarantee that important conformations are not overlooked. Therefore, unless a systematic method capable of always converging to the global minimum potential energy independent of the initial conformation is employed, the obtained minimum energy conformation will be limited by which initial conformations seemed appropriate to the researcher [3]. The need for a method that can GLOBAL MINIMUM ENERGY MOLECULAR CONFORMATIONS 139 guarantee convergence to the global minimum potential energy conformation is the motivation for this initial effort to introduce such a method for small molecules interacting with relatively simple force fields. 2. P r o b l e m D i s c u s s i o n The problem to be addressed in this work can be stated simply as follows: "Given the connectivity of the atoms in a molecule and the force field according to which they interact, find the molecular conformation(s) in the three-dimensional Euclidean space involving the global minimum total potential energy" The simplifications employed herein are as follows: 1. The molecular mechanical approximation of the Born-Oppenheimer surface is adopted. 2. Only palrwise, additive, two-body interaction terms are considered. 3. Covalent bond lengths and angles are assumed to remain at their equilibrium values. . The expression representing the pairwise non-bonded .atom potential interactions is assumed to be a function of only the Euclidean distance between the interacting atoms. Approximation (1) is well established in the study of molecular conformations providing tractable expressions for the total potential energy of the molecule. Simplification (2) is valid for nonpolar molecules whose most stable conformations are sought in this work. Assumption (3) is adopted for the sake of convenience and it is fairly accurate in most cases because covalent bond lengths and angles do not deform significantly from their equilibrium values without substantial increase in the potential energy of the molecule. Finally, assumption (4) is made routinely in the field of molecular mechanics without significant loss of accuracy. Under the aforementioned simplifications the expression for the total potential energy V of a molecule involves only the sum of a number of palrwise potential interaction terms. Each pairwise potential term is a function of only the Euclidean distances rlj between the interacting atoms which are directly related to the Cartesian coordinates xi, yi, zi of the atoms forming the molecule. Therefore, V can be fully represented in the coordinate space xi, y~, zi of the atomic coordinates. However, a number of equality constraints must be added in the formulation to reflect the fact that all covalent bond lengths and covalent bond angles are assumed to be fixed at their equilibrium values. After defining B to be the set of bonded atoms 140 M A R A N A S AND F L O U D A S and AfB the set of non-bonded atoms, V can be formulated in the Cartesian atomic coordinate space as follows: